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SELF CONSISTENT MODEL FOR THE EVOLUTION OF ECCENTRIC MASSIVE BLACK HOLE BINARIES 
IN STELLAR ENVIRONMENTS: IMPLICATIONS FOR GRAVITATIONAL WAVE OBSERVATIONS 

Alberto Sesana^ 
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ABSTRACT 

We construct evolutionary tracks for massive black hole binaries (MBHBs) embedded in a surround- 
ing distribution of stars. The dynamics of the binary is evolved by taking into account the erosion 
of the central stellar cusp bound to the massive black holes, the scattering of unbound stars feeding 
the binary loss cone, and the emission of gravitational waves (GWs). Stellar dynamics is treated in 
a hybrid fashion by coupling the results of numerical 3-body scattering experiments of bound and 
unbound stars to an analytical framework for the evolution of the stellar density distribution and 
for the efficiency of the binary loss cone refilling. Our main focus is on the behaviour of the binary 
eccentricity, in the attempt of addressing its importance in the merger process and its possible im- 
pact for GW detection with the planned Laser Interferometer Space Antenna {LISA), and ongoing 
and forthcoming pulsar timing array (PTA) campaigns. We produce a family of evolutionary tracks 
extensively sampling the relevant parameters of the system which are the binary mass, mass ratio and 
initial eccentricity, the slope of the stellar density distribution, its normalization and the efficiency of 
loss cone refilling. We find that, in general, stellar dynamics causes a dramatic increase of the MBHB 
eccentricity, especially for initially already mildly eccentric and/or unequal mass binaries. This affects 
the overall system dynamics; high eccentricities enhance the efficiency of GW emission, accelerating 
the final coalescence process. When applied to standard MBHB population models, our results pre- 
dict eccentricities in the ranges 10~^ — 0.2 and 0.03 — 0.3 for sources detectable by LISA and PTA 
respectively. Such figures may have a significant impact on the signal modelling, on source detection, 
and on the development of parameter estimation algorithms. 

Subject headings: black hole physics - methods: numerical - stellar dynamics - gravitational waves 



1. INTRODUCTION 

It is now widely recognized that massive black holes 
(MBH) are fundamental building blocks in the process of 
galaxy formation and evolution. MBHs are a ubiquitous 
compo nents of nearby galaxy nuclei (see, e.g., (Magorrian 
|et al. 1998), and their mas ses tightly corr elate with the 
properties of the host (Haring fc Rix [2004 and reference 
therein). In popular AcUM cosmologies, structure for 



namical simulations. However, whether viscous angular 
momentum extraction is ef ficient all the way d own to 
coalescence is questionable (Lodato et al. | |2009 ). 



mation proceeds in a hierarchical fashion (White fc Rees 
1978[ ) , through a sequence of merging events. It MBHs 



are common in galaxy centers at all epochs, as implied by 
the notion that galaxies harbor active nucle i for a short 
period of their lifetime ( Haehnelt fc Ree"s~| |l993), then 
a large number of massive black hole bmaries (MBHBs) 
are expected to form during cosmic history. The evolu- 
tion of such binarie s was firstly sketched by 'Be gelmaii^ 
B landlord fc Rees (1980), but after thirty years, sev- 
erai details ot the involved dynamical processes are still 
unclear. 

In stellar environments, the MBHB evolution proceeds 
via super-elastic scattering of surrounding stars in tersect 



In general, the vast majority of studies (mostly numer- 
ical) devoted to the subject have focused on the shrinking 
of the binary semimajor axis, because the relative small 
number of particles involved {N < 10^) make the eccen- 
tricity behaviour fairly noisy. However, eccentricity may 
play an important role in the final coalescence, because 
at a given semimajor axis, the coalescence timescale as- 
sociated to gravitational wave (G W) emission is much 



short er for very eccentric binaries ( Peters fc Mathews 
1963|. Moreover, having a trustworthy model tor the 



eccentricity evolution of the system may be of crucial 
importance for the practical detection of MBHBs in the 
forthcoming GW windows. 

MBHBs are infact expected to be the loudest sources of 
gravitation al radiation in the nHz -mHz frequency range 
(Haehnelt 199 4; JaflFc fc_Backer 2003; Wyithc fc Loeb| 
Enoki et al. 2004; Sesana et al 
112005; Rhook fc Wyithc 



ing the binary orbit (slingshot mechanism, Mikkola fc 
Valtonen~||1992 ), and the fate of the system depends on 



enet et al. 



2004, 200^ 
2005, Sesana. 



the supply ol stars available for such interaction. On 
the other hand, if the system is gas rich, torques exerted 
by a massive circumbinary disk have been prov en effi- 
cient in shrinking the binary dow n to '-^ 0.1 pc (Escala 



et al. 2005 Dotti et al 



2007), which is the current 



Vecchio fc Colaciiio||2008[ [Sesana, Vecchio fc Volonteri 
2009[ ) . The space-born e observatory Laser Interferometer 
Space Antenna (LISA, Danzmann et al."||1998 ) has been 
planned to cover the range of frequencies from 10~^ Hz 
to 0.1 Hz. Moving to the nanohertz f requency range, the 
Parkes Pulsar Timing Array (PPTA; |Manchestcr 2008) 



resolution limit ot dedicated smoothed particle hydrody 



^ Max Planck Institute for Gravitationalphysik (Albert Ein- 
stein Institute), Am Miihlenberg , 14476, Golm, Germany 



the Europ ean Pulsar Timing Array (EP'IA; 
al 



Janssen et 



] |2008[) and the North American Nanohertz Observa- 
tory to r Gravitational Waves (NANOGrav; Jenet et al. 
[] |2009[ ) are already collecting data and improving their 
sensitivity in the of 10^^ — 10^^ Hz window, and in the 
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next decade the planned Square Kilometer Array (SKA; 
Lazio [[2009 ) will provide a major leap in sensitivity. 

Besides tne technical progresses in the instrumenta- 
tion, the source signal modelling and the development of 
appropriate data analysis techniques for recovering the 
sources from the data stream are crucial for the success 
of the GW astronomy challenge. So far, most of the 
attention was focused on circular MBHBs. This seems 
to be justified, because GW emission is very efficient in 
dumpening the binary eccentricity, and since GW detec- 
tors (LISA in particular) are sensitive to the very end 
of the MBHB inspiral, sources are expected to be circu- 
lar when they enter the observable band. Consequently, 
most of the source modelling and the signal searches and 
analyses, e.g. the source inje ction in LISA mock data 
challenge ( Babak et al. |[2009 l, or the investigation car 
ried by the L ltiA parameter estimation task force ( Arun 



et. al 2009), relied on this assumption. 



However, both stellar and gas based shrinking mecha- 
nisms have proven to be efficient in increasing the binary 
eccentricity dQuinlan [19 96 Armitagc & Natarajan 2005; 
Sesana, Haardt fe Mada u 2006, 2008; Baumgardt, Gua- 
landris fc Porte gies Zwart 2 006; Ma tsubayashi, Makino 
& Ebisuzaki ' 2007 Bcrcntze n et 1 2009 ; Cuadr a et al. 
Miller k. Freitag 



[12009. ,Amar o-Seoane 
Seoane et al 



,2009, ^Amaro- 



2010 ), calling into question whether the 

assumption of circular orbits is justified for such GW 
sources. In this paper we construct a self-consistent sim- 
ple model for tracking the evolution of the MBHB ec- 
centricity (and semimajor axis) in stellar environments. 
We model the stellar distribution surrounding the bound 
binary as an isothermal sphere (p oc r~^) matching a 
cusp with a power law density profile p oc r^'^ inside the 
binary influence radius. The MBHB is evolved taking 
into account the scatt ering of bound stars leading t o the 
erosion of the cusp (Ses ana, Haardt fc Madau [2008 1, the 
subsequent scattering of unbound stars intersecting the 



Madau 



ni2006l), 

Mathew 



binary semim ajor axis dQuinlan ]T996| [Sesana, Haar dt fc 



and the efiicient GW emission stage (|Pe- 



lathews jl963 ) leading to final coalescence of the 



ters &: 

system. The main goal of the paper is to build sensible 
evolutionary tracks for the MBHB as a function of the bi- 
nary mass, mass ratio, initial eccentricity at pairing and 
cusp slope, and to show that viable MBHB evolution sce- 
narios predict a significant eccentricity in the frequency 
band relevant to GW observations with LISA and PTAs. 

The paper is organized as follows. In Section 2, we ex- 
tensively describe our model, defining the relevant phys- 
ical mechanisms and writing down the evolution equa- 
tions for the system. In Section 3 we provide further 
insights about the physics of the model, linking our treat- 
ment to the loss cone refilling theory. We present in de- 
tail our evolutionary tracks in Section 4, discussing the 
dependencies on the relevant model parameters, and we 
draw predictions for LISA and PTA observations in Sec- 
tion 5. Our main findings are summarized in Section 
6. 

2. INGREDIENTS OF THE MODEL 

2.1. Initial setup 

We consider a MBHB with mass M = M1+M2 {Mi > 
AI2) described by its semimajor axis a and eccentricity e. 
The system is embedded in a purely stellar background 



with a density profile described by a double power law, 
as follows: 



Pinf 



p{r) 



r > nnf. 



(1) 



Here ri^i is the influence radius of the binary, identifying 
the region where the gravitational potential is dominated 
by the two MBHs, and formally defined as the radius con- 
taining a stellar mass equal to M, and pinf is the stellar 
density at rinf . The density distribution is normalized to 
an isothermal sphere for r > rjnt, such a condition sets 



Pint 



2ttGi 



inf 



and 



nnf 



(3-7) 



GM 



0.8pc (3 



(2) 



(3) 



where Afg is the total mass of the binary in units of 
10^ McT), and we made use of the well established M — a 



relation in the form ( Tremaine et al. ]|2002 ) 



Me = 0.84a: 



70 



(4) 

((770 is the velocity dispersion in units of 70km s^^) to get 
rid of a in the last approximation. We identify the region 
r < Tinf as the inner cusp, and we use 7 = 1,1.5,2 cor- 
responding to nuclei characterized by cores/weak cusps, 
mild cusps and steep cusps respectively. 
N-body simulations of uneq ual mass binaries in the 



mass ratio range 0.01 — 10 (pB aumgar dt, Gualandris 
fc Port egies Zwart ] |2006[ [Matsubayashi, Makino &: 



Ebisuzaki 200T| have shown that dynamical friction is 



efficient in driving the secondary MBH much deeper than 
rinf in the potential well of the primary-cusp system. The 
two MBHs pair together forming a MBHB and continue 
to harden down to a separation at which the enclosed 
mass in the binary is of the order of M2 , without signifi- 
cantly affecting the stellar density profile. This is an indi- 
cation that the hardening is still driven by the dynamical 
friction exerted by the overall distribution of stars, rather 
than by close individual encounters with stars intersect- 
ing the binary orbit. In our model, we assume that M2 
is driven by dynamical friction down to a separation ao 
where the enclosed stellar mass in the binary is twice the 
mass of the secondary MBH 



ao = (3-7) 



GM 



q 



1 



1/(3-7) 



1 



q 



1/(3-7) 



(5) 

{q — M2/M1 is the mass ration of the binary system) 
without affecting the stellar distribution in the cusp sig- 
nificantly. At that point, three body interactions take 
over, and the binary evolution is dictated by individual 
encounters with stars intersecting its orbit. 

2.2. Physical mechanisms in operation 

On its way to final coalescence starting from oq, the 
binary is subject to three main dynamical mechanisms 
driving its evolution, namely: (i) the erosion of the cusp 
bound to the primary MBH, (ii) the scattering of un- 
bound stars supplied into the binary loss cone by re- 
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laxation processes once the stellar distribution is signifi- 
cantly modified by the MBHB, and (iii) the emission of 
GWs. The detailed description of each mechanism has 
been presented elsewhere, and the reader will be referred 
throughout this section to the appropriate references for 
further insights. The focus of the present work is to add 
them together coherently, to produce sensible, although 
admittedly very simplistic, evolutionary tracks for MB- 
HBs hardening in stellar environments. 

2.2.1. Bound cusp erosion 

The hardening of an unequal mass MBHB, with an 
initial semimajor axis ao and eccent ricity ep, in a bound 
cusp was extensively studied by Sesa na, Haardt fc Madau| 
|](p008), hereinafter SHM08. Using their formalism, two 
differential equations determine the rate of change of the 
orbital separation and eccentricity: 



da I 



2a2 



GM1M2 



da^dt 



de I 
di^'' 



A, 



d^N,; 

dat,dt 



da^, 



(6) 



(7) 



Here a* is the semimajor axis of a star bound to Mi 
and d^Nej/da^dt is the number of stars subject to sling- 
shot ejection in the semimajor axis and time intervals 
[a*, a* -I- da■^,], [t,t + dt]. The terms Ae, A£ are mea- 
sured from scattering experiments. The ejection rate 
d'^Ncj/da^.dt is instead computed by coupling the numer- 
ical results of the experiments to an analytic framework 
for the binary evolution, which is embedded in a cusp 
of the form described by equation Q, as detailed in 
SHM08. The major finding of SHMD8 is that the bi- 
nary hardens by a factor of ~ 10 by extracting the bind- 
ing energy of the stars in the cusp. During this process, 
e usually increases by a large factor, depending on the 
binary mass ratio and on the cusp slope. Results are 
tabulated in table 1 of SHM08. 

2.2.2. Slingshot of unbound stars 

The theory of MBHB hardening in a distribution of 
unbound field stars characterized by a density p and a 
velocity dispersion a was single d out by Quinlan ( 1996 1 
and ex tensively revisited by |Sesana, Haardt fe Madau 
(|2006|, hereinafter SHM06. The binary evolution can be 
expressed as a function of the dimensionless hardening 
rate H and eccentricity growth rate K as 



da I 
de I 



a^Gp 



aGp 



HK. 



(8) 



(9) 



The quantities H and K are related to the average en- 
ergy and angular momentum exchange between the stars 
and the binary in a single encounter, and are computed 
via extensive three body scattering experiments, as de- 
scribed, e.g., in SHM06. In general, hardening by scatter- 
ing of unbound stars becomes effective when t he binary 
reach the so called hardening radius, defined as (Quinlan 
[1|1996| 

ah « (10) 



This is the separation at which the specific binding en- 
ergy of the binary is of the order of the specific kinetic 
energy of the field stars. For a > ah, stars are basi- 
cally too fast to effectively exchange energy and angular 
momentum with the binary (soft binary regime); when 
a < ah, the binary tends to capture stars in short living 
weakly bound orbits, kicking them to infinity with v > a 
(hard binary regime). The transition soft/hard binary 
is rather smooth, and happens at about ah- Once the 
binary is hard, its hardening proceeds at about constant 
rate, as shown by the H tracks plotted in figure 3 of 
SHM06. Perfectly circular binaries tend to stay circular 
(because of the conservation of the Jacobian integral of 
motion in the 3-body problem), while even slightly eccen- 
tric binaries tend to increase their eccentricity, as shown 
by the K rates plotted in figure 4 of SHM06. 

2.2.3. Gravitational wave emission 

For our purposes, the effect of GW emission can be 
modelled in the quadrupole approximation. Under this 
assumpt ion, the evolution equation s for the system are 
given by Peters & Mathews ( 1963 1 



da I 
dt 's^ 



de I 
dt 



64 M1M2M 
'y^a3(l-e2)7/2 



1 



73 

24^ 



37 

96^ 



64 G^ M1M2M 



Fie) 



304 M1M2M 



15 c' 



a4(l _ e2)5/2 



1 



121 
304' 



(11) 
(12) 



The fun ction F{e) is defined by the last equality in equa- 
tion 11 The shrinking rate is a strong factor of a, mean- 
ing that GW-driven hardening is effective only at small 
separations. The eccentricity evolution rate is also a 
strong function of a and e itself, and it is always negative. 
GW emission, therefore, is very effective in circularizing 
MBHBs, which, in turn, is the reason why little attention 
has been paid so far to eccentric systems in the context 
of GW detection. 

2.3. General equations for the binary evolution 

Having identified the relevant mechanisms at play, we 
can put the pieces together by writing the evolution of 
the binary as 



dt 

de 
dt 



dt 

Ede I 



(14) 



where i = b,u, gw labels the three mechanisms consid- 
ered. Here we handle MBHBs in stellar environments, 
and the relevant scale of the stellar distribution is de- 
fined by the sphere of influence rinf of the massive black 
hole binaries. It is then natural to consider as relevant 
parameters the stellar density and velocity dispersion at 
the influence radius, p-mf and Uinf — a (for an isothermal 
sphere the velocity dispersion is independent on radius). 
It is then instructive to recast the MBHB dynamics in the 
dimensionless H and K formalism proposed by Quinlan, 
to compare the dimensionless rates given by each mech- 
anism. The global evolution of the system can be then 
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written as a generalization of equations ([s]) and (|9| in 
the form 



da 
~di 

de 
dt 



inf 



aGp 



inf \ ^ 



where Hi is triviaUy defined as 

a da I 



a'^Gpinf dt 



(15) 



(16) 



(17) 



and 



K, = a- 



de I 
dt ' 



E 



da I 



We integrate the coupled difi'erential equation s (jlSj ) 
( 16 ) starting from ap. As described in Section 27l (e 



(18) 

and 
equa- 
tion ([5|) ) , the value of oq is set by the total mass of the 
binary M, the mass ratio q = M2/M1, the cusp slope 
7 and the stellar velocity dispersion a. In our default 
models we force a to obey the M — a relation given 
by equation Q. In this manner there is a one-to-one 
correspondence between the mass of the binary and a, 
i.e., equally massive binaries are embedded in identical 
isothermal spheres. Having set the normalization of the 
stellar distribution and the initial separation oq, the evo- 
lution of the binary depends on the four parameters Mi , 
q, Co and 7. We extensively sample this parameter space 
as following: 

. log(Mi/M0) = 2, 3, 4, 5, 6, 7, 8, 9, 10, 11 

• g = 1,1/3,1/9, ...,1/729 

• eo = 0.01,0.1,0.3,0.6,0.9 

• 7 = 1,1-5,2 



for a grand total ofl0 x7x5x3 = 1050 simulations. 
Even though the binary evolution under the effect of stel- 
lar encounters is basically scale free, the simulation of 
systems with different absolute masses was necessary to 
match together the scattering-driven phase to the GW- 
driven phase, which instead is highly mass dependent. 
The assumption of different eccentricities at the moment 
of pairing takes the environmental effects affecting the 
dynamical friction stage into account. In general, during 
the merging process, galaxies capture each other on a 
very eccentric orbit, which is reflected in the initial tra- 
jectories of the two MBHs (still at kpc separations at this 
stage). Dynamical friction against massive, rotationally 
supported, circum binary disks has been proven to cir- 



cularize the orbit (Dotti, Colpi & Haardt 2006). How- 



ever, this is not in general true in gas poor environments, 
where the pr ocess is driven by interaction with t he stellar 
distribution ( Colpi, Mayer &: Governato~||1999 1, and the 
eccentricity of the MBHH at the moment of pairing may 
retain memory of its initial value, or may, in general, be 
different than zero. 

We also consider four alternative models to address the 
impact of the assumed M — a relation and the choice of 
normalizing the efficiency of unbound scatterings to pinf . 



Let us denote with a the value of the velocity disper- 
sion predicted by the M — a relation for a given M. We 
consider models with velocity dispersions equal to 0.70" 
and 1.3fT, which is approximately the range of variance 
of a fo r a given MBH mass measured in the M — ct re- 



lation (Haring & Rix 



2004j ). We also run models with 
two ditrerent normalizations for the unbound scattering 
process: a fast model normalized to lOpinf and a slow 
model normalized to O.lpinf. The mot ivation for this set 
of runs will be clarified in Section 13.21 

To practically evolve the binary we make use of the 
results of the scattering experiments with unbound and 
bound stars performed in SHM06 and SHM08. In those 
papers, the quantities Ae, A5, d^Ng-Jda^dt (for the 
bound scatterings), H and K (for the unbound scatter- 
ings) were recorded on a grid of a and e, covering the 
relevant dynamical ra nge . Th e ev olution of the coupled 
differential equations ( 15 1 and ( 16 ) is performed by inter- 



polation over the grid as the binary evolves. In SHM06, 
unbound scatterings were carried out for all the consid- 
ered mass ratios down to q = 1/243. To complete the 
sample we carried out additional experiments for the case 
q = 1/729. SHM08, instead, focused on unequal MBHBs, 
with g<0.1. To complete the mass ratio sample, we ran 
experiments for the cases q = 1, 1/3. We then have all 
the bound and unbound scattering experiments results 
spanning the q and bq range of interest. 

2.4. Limitations and caveats 

Our evolutionary tracks are computed in a self- 
consistent way, summing together the effects of differ- 
ent mechanisms. However, we should be aware of the 
several limitations and simplifications we have adopted. 
One major caveat is the extension of the bound scatter- 
ing experiments to mass ratios of g > 1/9. It is infact 
unlikely that, in such cases, M2 would reach without 
affecting the stellar cusp at all. Cusp disruption would 
start earlier (especially in the equal mass case), and the 
cusp erosion phase described here may not be a trustwor- 
thy description of reality. However, we find that the im- 
pact of the bound cusp erosion on the binary evolution is 
smaller for larger mass ratios. This is because for g — >■ 1, 
oo ^ 7'inf , (see figure [1]) and the impact of the binding en- 
ergy extraction in the total energy budget of the system 
becomes less significant. The eccentricity evolution in the 
cusp erosion phase is only mild when q — 1,1/3, implying 
that our approximate treatment would not significantly 
affect the overall results. Another caveat to bear in mind 
is that the three body scattering is a scale free problem 
as long as ^ M2 ■ Even though we compute 'a poste- 
riori' evolutionary tracks for systems with Mi = 100 M0 
and q = 1/729, we consider our results meaningful only 
when M2 > 100 M0. This is also reasonable, since we 
are interested in the evolution of MBHBs. Moreover, if 
M2 is small (< 10* M©), the amount of stars interacting 
with the binary is also quite small, i.e. the granularity of 
the problem increases. Our smooth evolutionary tracks 
should then be interpreted more as 'trends', or 'mean 
evolutions' rather than paths followed by each individ- 
ual binary. We have also not included the possibility of 
stellar tidal disruption, which has been shown to be an 
efficient process in the cusp erosion phase, especially in 



the mass ratio range 0.01 < q < 0.1 (Chen et al. |2009 ) 
In general, the inclusion of tidal disruptions mildly en 
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FlG. 1. — Relevant lengthscales in the binary hardening problem 
as a function of the binary mass ratio for two different absolute 
values of the binary mass M = 10^ M0 (thick lines) and M = 
10^ M0 (thin lines). For each set of curves, from left to right, 
we plot rinf,ao,<ih and a™ (for a circular binary, i.e. assuming 
F(e) = 1, see equation | |11[ |), as labelled in the figure. The cusp 
slope is fixed to 7 = 1.5. 



2010[ ), be- 

ly have a* < a, i.e. they 



hance the eccentricity increase (|Chen et al. 
cause disrupted stars preferentia 

would drive binary toward circularization if ejected (see 
SHM08 for a detailed discussion of this effect). Lastly, 
there is a somewhat net distinction between bound and 
unbound scatterings in our formalism, which is certainly 
oversimplistic, since in reality relaxation processes will 
mix-up the different stellar populations. The appear- 
ance of distinctive features in the transition between the 
bound and the unbound regime can be therefore consid- 
ered somewhat artificial, the reality would probably be 
more gentle. We do not believe that this has a major 
impact on our main results. 

3. PHYSICS OF THE BINARY EVOLUTION 

The relative contributions of the three mechanisms 
considered in the previous section are set by the typi- 
cal lengthscales below which they become effective. As 
stated before, cusp erosion becomes effective at oq. At 
this point iJf, and Ki, start to dominate the MBHB evo- 
lution. 

On the other hand, scattering of unbound stars be- 
comes fully effective at the hardening radius Uh defined 
by equation ( 19 ) 



1 



ah 



q 



4(3-7) \i + q. 

q 



(2-7)/(3-7) 



ao 



^0.2pcAfg/^ 



1 



(19) 



We see that — (l/4)ao independently of the mass ra- 
tio for 7 = 2, and in general, with decreasing 7, the ratio 
tth/ao becomes smaller and q dependent. As the binary 
shrinks, the central cusp is depleted and the scattering 
of unbound stars refilling the binary loss cone (i.e., the 
family of stellar orbits intersecting the binary semimajor 



axis, see next section) becomes dominant. In this sec- 
ond stage, occurring at approximately O.lao the binary 
evolution is determined by and K^- 

Gravitational radiation will eventually take over at 
Ogw, driving the binary to the final coalescence. 

Q>gw can 

be derived by imposing da/dt\i^ — da/dt\g^ ; rearranging 
and dividing by au we find 



1287r(3 - 7)F(e) 



1/5 



0.00027pc (3 - -if'^Fiefl^Ml'^ 



^1/5 



(l+g)2/5 



(20) 



where F{e) is defined by equation (11), and we again 
used the M — a relation to write the last approximation. 
Figure [T] highlights the behaviour of the lengthscales of 
the system (nnf , ag, a/j, agw) as a function of q for two 
selected values of M and 7 = 1.5. It is easy to see that, 
in general, Cgw <C ah- For example, assuming circular 
binaries with q — I, H — 15, and a = lOOkm/s (or 
M = 4x 10^ Mq, according to the M — a) gives a^-^/ah ~ 
2x10^''. flgw can approach ah if the mass ratio is extreme 
[q ^ 10~^) and the velocity dispersion is very large (e.g. 
if the binary is massive; a > 300km/s or M > 10^ M©), 
as shown by the set of thin lines in figure [l] In general, 
figure [1] clearly shows that there are three well defined 
zones where one single shrinking mechanism is dominant 
among the others, and we have: 

flgw < ah < ao < nnf- (21) 

3.1. H and K formalism in the loss cone framework 

It is instructive at this point to frame the hardening 
rate given by equation ([sl in the context of loss cone re 
filling theory ([F rank & Rccs 1976 Lightman & Shapiro 



1|1977 |Cohn fe'Kulsrud J978f ~Mter the binary has de 



pleted all the stars intersectmg its orbit (formally defin- 
ing what is called the 'binary loss cone' of the s tellar dis- 
tribution function, see Milosavljevic fc Merritt | (2003) for 
a comprehensive review), its further hardening depends 
on the rate F at which such orbits (i.e the loss cone) are 
refilled. Loss cone refilling can proceed by diffusion of 
stars either in energy (e*) or in angular momentum (j*). 
In general, stellar encounters are much more efficient in 
changing the star angular momentum, and diffusion in 
the j* space is the relevant process. The loss cone refill- 
ing depends on the average Aj* experienced by a star on 
a almost radial orbit during one orbital period. If this 
change is larger than the size of the binary loss cone in 
the angular momentum space, j\c ^ V2GMa, then stars 
are easily scattered back and forth into the loss cone, and 
the loss cone is filled. In this regime, the supply rate of 
stars from a given distance to the binary r is given by 
( [Lightman fc Shapiro ||1977| [Perets fc Alexande7]|2008[ ) 



dVf _ aN,{< r) 



dlogr 



P{r) 



(22) 



Here P{r) is the typical period of a star on an almost 
radial orbit coming from a distance r and Ni.(< r) is 
the number of stars enclosed in a sphere of radius r 
around the MBHB. Let us focus on stars coming from 
r > Hnf- Assuming an isothermal sphere iV*(< r) = 
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(M/m*)(r/rinf ), and that the typical period of a sta r on 
a radial orbit is P{r) ^ r/a, integrating equation (22 1 
from Tinf to c», and using equation ([2]), we get 

r,«2.^^. (23) 
TO* a 

Let us contrast this result with equation In the 
Quinlan formulation, the binary is embedded m a homo- 
geneous stellar field with density p. The hardening rate 
H is then derived writing the interaction rate as a flux 
of stars through the binary cross section, namely 



r 



Q 



(24) 



where p/m^, is the number density of field stars, v their 
velocity at infinity (i.e., far from the binary) and S the 
binary cross section. If b is the star impact parameter at 
infinity, and if we assume the encounter to be relevant 
only for b < 6max, then S = Trb^^^. Relating b to the 
maximum approach x to the binary via gravitational fo- 
cusing (6^ — 2GMx/v'^), and replacing Xmax = o, (stars 
have to cross the binary semimajor axis to exchange en- 
ergy and angular momentum efficiently), we get 



To = 2n 



M Gpa 



(25) 



By comparing equation (23 1 and (25), if we identify the 
intruder velocity at infinity v with the dispersion velocity 
in the isothermal sphere cr, we see that our Hu and 
prescriptions for the scattering of unbound stars corre- 
spond to considering the loss cone always full at rinf. 

3.2. Loss cone refilling 

By normalizing Hu and Ku to a and pinf, our model 
implicitly assumes that the loss cone is always full at 
''' > Hnf (i-e., in the so called pinhole regime), and 
empty for r < rinf (i.e., i n the so called diffusive regime, 

The status of the loss cone 

then enter as a parameter in our formulation, set by 
our choice of normalizing the system to pi^f. The is- 
sue of what is the physical mechanism that keeps the 
loss cone full at r > rini is not addressed in this paper. 
We will only briefiy discuss here the plausibility of such 
scenario. After the loss cone is depleted, in absence of 



any other physical mech anism, two body relaxation (Bin- 



ney &: Tremaine 

refilling 

in real galaxies (Merritt & Szell 2006), and under this 



1987 ) sets the timescale for loss cone 

i'his I S usually longer than t he Hubble times 



assumption, the loss cone is, in general, in the diffusive 
regime way beyond rjnf. However, in more realistic situ- 
ations, a myriad of other physical factors play a substan- 
tial role, shortening the loss cone refilling timescale. It 
has been shown that axisymmetry and in particular tri- 
axialit y ( |Yu ||2002yMerritt k Poon ||2004||Berczik et al. 
|2006 1 are very ettective in repopulatmg the loss cone, 
driving stars on chaotic and centrophilic orbits toward 
the black hole. Moreover, the matter distribution in stel- 
lar bulges is far from being smooth. Inhomogeneous con- 
centrations of matter, such as massive star clusters or gi- 
ant molecular clouds (the so called massive perturbers), 
have been proven efficient in perturbing stellar orbits, 
significantly sh ortening the loss cone refilling timescale 
( Perets fc Alexander ,2008| ). We should note that these 



mechanisms are likely to operate efficiently for r > rjnf. 
Inside the MBH influence radius, indeed, the presence 
of a central massive object dominating the gravitational 
potential tends to force the stellar system toward a more 
spherical symmetry (and the triaxial structure may be 
erase d in the central region, see e.g [Merritt fc Quin~ 



lanl (1998)). For the same reason, massive perturbers 
tend to be stripped by the MBH tidal field, and would 
hardly survive in this region. Also the two body relax- 
ation timescale increases for decreasing r if the potential 
is dominated by a central object, because the velocity 
dispersion of the system, which in this region is the Ke- 
plerian velocity around the MBH, scales as r~^/^, and 
the relaxation times cale has a strong dependenc e on the 



1987[ ). More- 



velocity dispersion (Binney & Tremaine 
over, inside rjnf, the ejection of stars bound in the cusp 
depletes those orbits with energy e* < GM/(2rinf), and 
since e* diffusion is much less efficient than j* diffusion, 
the contribution to the loss cone refilling coming from 
T < ''inf should be, in general, negligible. It is then rea- 
sonable to assume a full loss cone for r > rjnf, and other- 
wise empty. We shall test the implication of such hypoth- 
esis by v arying the normalization used in equations (15) 
and (16). We therefore test models normalized to lOpinf 
and 0.1 Pint, which correspond to consider the loss cone 
full only for r > lO^/^nnf , or down to r = 10^('^^''''/^7'inf 
(where 7 is the slope of the inner cusp as defined by equa- 
tion ([!])) respectively. The consideration of smaller values 
of the relevant density assumed for the diffusion process 
(O.lpinf) serves also as a test for the robustness of our 
results against different outer density profiles. Although 
our refilling rate is derived for an isothermal distribu- 
tion, that reasonably fit the measured density profiles of 
severa nearby galaxies and of the M ilky Way (Lauer et 



al. 



les s 



1995 Dehnen & Binney 1998), ma ny other galax- 



low shallower outer density profiles (Ferrarese et al. 



1 |1994[ [Gebhardt et al ||1996| [Rest et al. ||200l[ ) and 



the bulk of stars participatmg to the refilling mechanism 
may come from slightly larger radii, where the density is 
lower. We will show (section 5.2 and figure 7) that this 
have a minor impact on our results, changing the eccen- 
tricity of the system in the observable GW bands by a 
factor <2. 



4. EVOLUTIONARY TRACKS FOR MBHBS 

In this section we present the MBHB evolutionary 
tracks produced by our hybrid model. As discussed in 
the previous section, the binary goes through three sub- 
sequent phases, which are in general distinct. In the 
discussion we will simply identify them as the bound 
phase (erosion of the bound cusp), unbound phase (scat- 
tering of unbound stars refilling the loss cone) , and GW 
phase (where GW emission become more efficient then 
the unbound scattering) . Each phase is characterized by 
its proper Hi and Ki. In the following compilation of 
plots we will present the evolution of the global rates 
H = J2i Hi and K = Ki, it will be clear by looking 
at the figures which particular mechanism dominates in 
each region of the binary evolution. Each of the figures 
[2] [3] and [4] shows the quantities e(t), a/ao{t), K{a/ao) 
ana _ff (a/oo) for different q, as labeled in each panel. In 
the discussion we will simply refer to the panels as e, a, 
H and K panels. 



Eccentric massive black hole binaries 



7 



10* 10= 10* 10' 10= 10* 10^ 10* 10' 10= 

= I iiiii^ I iiiiiiii I iiiiiii I iiiiiiii y 1 



1 o.io.oiio-no-+io-» 1 0.1 o.oiio-»io-''io-* 



0,2 Hiiii I piiii I piiiii I piiiii I piiii 

0.1 I- q=i -i 

— ■ ^ 



0,4 
0.2 





10-3 r q=V8i \i _; 



10" 10= 10= 10' 10= 10* IflS 10* 10' 10= 
t[yr] 



1 0.1 0.01 10-3 lO-* 1 



0.1 0.01 10-3 10-4 



Fig. 2. — MBHB evolutionary tracks produced by our model by assuming 7 = 1.5 and eq = 0.1. Left plot: in each of the four pairs of 
panels we plot the eccentricity (top) and semimajor axis (bottom) evolution as a function of time. Different panels refer to different values 
of q, as marked by the inset labels. Different linestyles correspond to different binary masses: M = 10^ Mq (solid), 10^ Mq (dotted), 
10^ M0 (short— dashed), 10* Mq (long-dashed), 10^ M0 (dotted-dashed). Right plot: corresponding binary eccentricity growth rate (top 
panels in each of the four sectors) and hardening rate H (bottom panels in each of the four sectors), as a factor of the binary separation 
normalized to ag. Linestyle as in the left plot. 
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Fig. 3. — Same as figures but now assuming M = 10^ Mq and varying the initial eccentricity bq. In each individual panel different 
linestyles are for eo = 0.01 (solid), 0.1 (dotted), 0.3 (short-dashed), 0.6 (long-dashed), 0.9 (dotted-dashed). The inner cusp slope is fixed 
to 7 = 1.5. 



4.1. Dependence on the model parameters 

The evolution of the system depends on the chosen 
values for the parameters Mi, q, eo and 7. Such depen- 
dencies are extensively illustrated in figures |2j [3] and |4j 
Let us consider each parameter separately, by starting 
with those defining the masses of the system: Mi and q. 

The evolution of the MBHB as a function of Mi and 
for different q is plotted in figure [2j assuming 7 = 1.5 and 



Ci — 0.1. Since our treatment of stellar scattering is scale 
free, the value of Mi affects the system evolution only, by 
setting the relative gap between and Ogw, which has a 
mild dependence o n M -\ . By substituting the M — a re- 
lation in equation (20l we have infact a^/ag^ oc 
i.e. the gap is larger for lighter systems. This means 
that, in general, lighter binaries become more eccentric, 
because, after the bound scattering phase (which is ba- 
sically scale free by construction) they evolve under the 
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Fig. 4. — Same as figurcl2|but now assuming M = 10® Mq, eg = 0.1 and varying the cusp slope. Different linestyles are for 7 = 1 (solid), 
1.5 (dashed), 2 (dotted— dasnod) . 



effect of unbound scattering for a larger portion of their 
dynamical range, as it becomes clear by looking at the K 
panels of figure [2] For equal mass binaries with eg = 0.1, 
the eccentricity grows only to 0.15 for Mi = 10^ M0, 
and up to 0.3 for M\ = 10^ Mq. The mass dependence 
of the eccentricity growth is also evident for binaries with 
q = 1/9, while it tends to disappear for lower mass ra- 
tios where the eccentricity evolution is dominated by the 
bound phase. The q dependence of the eccentricity evo- 
lution is emphasized by the four different quadrants of 
figures |2j |3] and [4] Let us consider again figure [2j The 
main result here is that the eccentricity growth in the 
bound phase, at least when bq is small, is in general much 
larger for lower values of q, as explained in detail in Sec- 
tion 4.1 of SHM08. This is nicely shown by the K panels: 
as q decreases from 1 to 1/729, the peak in the K rate in- 
creases from ^ 0.05 to ~ 0.6. The value of q also sets the 
relative weight of the bound and of the unbound phases 
in the hardening process. Although a^/ah is only mildly 
depe nde nt on q (depending on the cusp slope 7, see equa- 
tion ( 19 )), tth/cLgw is a strong function of q, (see equation 
(20)) and it can be even less then one for high Mi and 
small q (as shown in figurejl]). Therefore, as q decreases, 
the eccentricity growth is dominated by the bound phase. 
For eo = 0.1, the maximum eccentricity reached by the 
binary at the end of the stellar driven phase increases 
from - 0.2 for g = 1 to - 0.9 for q = 1/729. The trends 
with Ml and q presented in figure [2] are preserved when 
changing 7 and cq. 

The dependence of the MBHB evolution on cq is stud- 
ied in figure [3] where we fixed Mi = 10^ Mq and 7 = 1.5. 
The binary eccentricity, in general, tends to increase in 
the scattering phase regardless on the value of eg. When 
Co = 0.01, binaries with q > 0.1 experience only a mild 
increase in their eccentricity, up to a value ~ 0.2, while 
binaries with smaller q can reach e > 0.8. Binaries with 
Co > 0.3 tend to reach eccentricities larger than 0.9 re- 
gardless on q, Ml and 7. The evolution of H is basically 



unaffected by the eccentricity of the system in the bound 
and unbound phases (in general the average energy sub- 
tracted to the binary by the star is not affected by the 
binary eccentricity, see e.g. SHM06), while K is interest- 
ingly larger for higher eo when q is large, and viceversa, 
it decreases with increasing eg for small q. 

The impact of the assumed slope 7 of the density pro- 
file is highlighted in figure |4] for binaries with Mi = 
10® Mq and Cq = 0.1. In general, MBHBs in steeper 
cusps evolve faster, but to a lower maximum eccentricity, 
during the scattering phase. As explained in SHM08 this 
is because, in the scattering process, stars with a* > a 
tend to increase e, while stars with a* < a tend to de- 
crease it, and the relative weight of the formers is larger 
in shallower cusps. Moreover, by increasing 7, the dy- 
namical range covered by the scattering process is much 
shorter (especially for low g), because the ag — Ogw gap is 
smaller, and there is less room for significant eccentric- 
ity growth. The K rates are mildly affected by 7, with 
higher value of 7 resulting in smaller K, as explained 
before. Also notice that the absolute value of H in the 
bound phase increases a lot with 7. This is because the 
value of Go is much smaller for high 7, and the shrink in 
the bound phase is accordingly much faster. In general, 
for any value of 7, the eccentricity reached at the end 
of the star scattering phase is > 0.7 for q < 0.1, while, 
again, equal mass binaries experience a less pronounced 
increase in the eccentricity. 

Figures [2] [3] and |4] allow also a detailed study of the 
evolutionary timescale as a function of the system pa- 
rameters. Since the bound phase is usually much faster 
than the unbound one, the evolution timescale of the sys- 
tem is set by a/d — a/{GpHa) oc 1/a. The coalescence 
timescale is then set by the bound-GW transition occur- 
ring at agw By substituting Ogw given by equation (20 1 
in the timescale definition above, we get for the coales- 
cence timescale 

Tc OC F(e)~i/5(3 - 7)^/^(7^1/^(1 + g)2/^ (26) 
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Fig. 5. — Contour plots of e computed at Jlisa = 5 x 10 ^Hz, in the {Mi,q) plane for selected model parameters, as labelled in each 
panel. Note that we excluded the bottom-left region, corresponding to systems for which M2 < 100 Mq . 



Firstly, we notice that Tc is independent on the absolute 
value of the MBHB mass, in agreement with figure [2j 
This is a consequence of normalizing the stellar distri- 
bution outside Tinf to an isothermal sphere obeying the 
M — fj relation. Tc has also a very mild dependence on q, 
increasing by a factor of ^ 3 when the mass ratio drops 
from q — \ioq~ 1/729, as shown by the correspondent 
panels in figure [2] High values of the maximum eccentric- 
ity accelerate the coalescence by a factor i^(e)^/^, which 
is ~ 5 for e = 0.9, this effect is clear in the a and e panels 
of figure [3] The impact of 7 is also quite mild, in spite 
of the 9/5 exponent, and it modifies Tc by a factor of 
~ 3 — 4 (see a and e panels in figure |4]). In general, we 
find lO^yr < < few x lO^yr. 

4.2. Comparison with numerical works 

The evolution of MBHBs in stellar environments has 
been tackled by several authors by means of full N-body 
simulations ( Milosavljevic fc Merritt |2001 Hemsendorf 



2006 


Matsubayashi, Makino & Ebisuzaki 2007 


Mer- 


ritt, M 


[ikkola & Szell 1 


20071 IBerentzen et al. 


120091 


Amaro-Seoane, Miller & Freitag | 


20091 lAmaro-Seoanel 


et al. 2010|). However, the limited number of particles 



Sigurdsson fc Spurzem ||2002[ [Aarsetb ||2U03| |M akino fe' 
Funato .2004;,Baumgardt7C!ualandris fc Portegies Zwart, 



haviour for the binary eccentricity, and it is difhcult to 
draw conclusions about the general trends behind the nu- 
merical noise. We can compare our results with N-body 
simulations carried out in two regimes: q — 1 (equal mass 
inspiral s) and q = 1/1000 (interinediate MBH-MBH in- 
spiral). Milosavljevic fc Merritt | (2001) carried out nu- 
merical integration of equal MBHBs embedded in two 
merging isothermal cusps (7 = 2). Starting with cir- 
cular orbits they find a mild eccentricity increase to a 
value of <0.2 during the stell ar driven hardening phase, 
consistent with our findings. [Merritt, Mikkola fc Szel^ 
^1 (|2007| considered equal MBEIBs embedded in Dehnen 
density profiles ( [Dehnen ^1993) with 7 = 1.2 with differ- 
ent initial eccentricities. Again, they find that circular 
binaries tend to stay circular, while eccentric binaries 
tend to increase their eccentricities in reasonable agree- 
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Fig. 6. — Same as figure [s] but now for /pxA = 5 X 10 ^Hz. Note the different scale in the x-axis, dictated by the fact that PTAs are 
sensitive to more massive systems emitting at lower frequencies. 



mcnt with the prediction of scattering experiments, and 
consequently, with the tracks we presented i n figure [3 
for th e q — 1 case. Simulations carried by [Aarseth 
(|2003|) and Hamsendorf et al. (2002) produce MBHBs 
0.8 at the moment of pairing, with e sub- 
again consistent with 
studied the evolu- 



in the next section. 



by Amaro-Seoane et al 



Sim i lar co nclusions are reported 
( |20T(| 



witTi eg 

sequently increasing up to >0.95 

p09p 



On the small q side. 



sini ulatio ns were performed by iBaumgardt, Gualandris 
fc Portcg ies Zwart | ( 2006 1 and [Matsubayashi, Makino 
& Ebisuzaki 
T77^ 



our findings. Berentzen et al. 

tion of equal MBHBs in rotating systems described by a 
King stellar distribution. They also find quite eccentric 
binaries at the moment of pairing (e > 0.4), and the sub- 
sequent evolution leads to eccentricities larger than 0.95 
at which point GW emission takes over, again consis- 
tent with our findings. [Amaro-Seoane, Miller &: Freitag| 
d ( |2009| ) focused on intermediate MBHBs (M ~ lO'^M©) 
m massive star clusters. They employ a machinery sim- 
ilar to ours, coupling full N-body simulations to three 
body scattering experiments. Their binaries have sig- 
nificant eccentricity 0.5 — 0.6) at the moment of 
pairing, and the predicted range in the LISA band is 
0.1 < e < 0.3, in good agreement with the results shown 
by the eccentricity maps in figure [s] that we will describe 



assuming a stellar density profile 
7 = i.YO. When properly rescaled, the eccentricity in- 
crease found in both papers agrees surprisingly well with 
our predictions based on the hybrid cusp erosion model. 
Unfortunately, we did not find any mention in the N- 
body literature about the eccentricity evolution for in- 
termediate values of q (0.01 < q < 0.1), and it would be 
useful to compare our results with N-body simulations 
in this intermediate range. We find, however, the overall 
good agreement, at least in the trends, shown at the two 
extremes of the mass ratio range comforting. 

5. ECCENTRICITY IN THE LISA AND PTA WINDOWS 

One of the main goals of the present study is to draw 
sensible predictions for the eccentricity of MBHBs emit- 
ting GWs in the LISA and in the PTA frequency ranges. 
So far, most of the work related to source modelling. 
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signal analysis and parameter estimation relied on the 
assumption of circular orbits. This seems reasonable be- 
cause GW emission is very efficient in dumpening the 
binary eccentricity, and since GW detectors {LISA in 
particular) are sensitive to the very end of the MBHB 
inspiral, sources are assumed to be circular when they 
enter the observable band. However the level of residual 
eccentricity critically depends on how large e is at the 
transition between the stellar hardening and the GW 
phases. In the scenario proposed here, such value can 
easily be larger than 0.9, implying non negligible resid- 
ual eccentricities in the frequency bands to be probed by 
future GW detectors. 

5.1. Eccentricity maps 

To convert our evolutionary track into predictions for 
GW observations, we proceed as follows. For each 
MBHB (uniquely defined by Mi, q, 7, cq), we convert 
the a/ao(t) tracks into a(t) tracks, and then we com- 
pute fk{t), the orbital frequency of the binary, simply by 
assuming Kepler's law. Having e{t) and fk{t) we then 
construct the e(/fe) evolution from the moment of pair- 
ing to the final coalescence. We then select two frequen- 
cies appropriate for LISA and PTA campaigns, Jlisa 
and /pTA respectively, and evaluate ej/^j-^^ and e|/p^^. 
Remember that for circular binaries, gravitational radia- 
tion is emitted at /gw = 2/^. We pick Jlisa = 5 x 10"^ 
Hz (corresponding to /gw = 10"'' Hz, which is approxi- 
mately the lower bound of the LISA band) ; on the other 
hand, we use /pta = 5 x 10"^ Hz (corresponding to 



10 Hz, which is approximately the frequency 



at which a 5-to-lO yr PTA campaign will be most sen- 
sitive to). The results are shown in figures [5] and [6] as 
contour plots e\f^jgy^{Mi, q) and e\fp.j.f^{Mi,q), for se- 
lected value of 7 and Bq, as labelled in the figures. Let 
us start discussing the LISA case. Firstly, we limited Mi 
to an upper value of 10^ Mq, since the inspiral of bina- 
ries with higher masses will fall outside the LISA band 

I We also excluded from our contour plots systems with 
2 < 100 Mq , for the reasons discussed in Section |2.4 



The general trend is that lighter unequal mass binaries 
tend to have larger e when they enter the relevant fre- 
quency range. The mass trend is easily explained by the 
fact that our treatment is largely mass invariant, but the 
absolute frequency of the system is not! Same stages of 
the MBHB evolution correspond to progressively lower 
frequencies as Mi increases, and since we are in the GW 
dominated phase (and thus de/df < 0), more massive 
systems have lower eccentricities at a given frequency. 
Binaries with smaller q are in general more eccentric 
because the eccentricity growth they experience in the 
stellar scattering phase is larger. If binaries are approx- 
imately circular at the moment of pairing (cq — 0.01), 
then the maximum eccentricity in the LISA band is ~ 0.2 
when Ml ^ 10"* M© and q<0-l, and the general trend 
is largely independent on 7. This is the result of two 
competitive effects: milder cusps lead to larger values of 

^ The higher frequency signal coming from the coalescence and 
ringdo wn as w ell as h igher harm onic corrections to the late inspiral 
phase i jPorter fc Cornish ||2008t are likely to push the detectable 
mass limit close to 10° Mq; however, the imprint of any residual 
eccentricity would be very small and hard to observe for such ex- 
treme masses. 



e, but in this case GW takes over earlier (because the 
timescale of the 3-body scattering evolution is set by the 
density at rinf which is lower for milder central cusps, be- 
ing Tinf itself larger) , and it has more time to circularize 
the orbit before the binary get to Jlisa- If binaries are 
significantly eccentric at the moment of pairing (eg = 0.6 
as a study case), then the q dependence in the contour 
plots almost disappear, because the eccentricity growth 
in the scattering phase is very efficient irrespective of q 
when eo is large. In this case, light MBHBs (Mi < 10"*) 
may reach Jlisa with eccentricities up to ~ 0.5. 

The situation is even more 'dramatic' for PTA ob- 
servations. PTAs are sensitive to much larger masses 
{Ml > 10^ M0) emitting at much lower frequencies 
(/ ~ 10~^ Hz). Systems are, in general, caught far from 
coalescence (the typical t ime to coalescence is ~ 10"* yr, 
Sesana & Vecchio 2010) and they did not have much 
time to circularize under the effect of GW emission. Be- 
cause of this, even if binaries were circular at the moment 
of pairing, eccentricities can be as high as 0.7 in the PTA 
band, with the same trend observed for the LISA case 
(i.e., lighter binaries with smaller q are more eccentric). 
If Co = 0.6, then all the systems with Mi < 10^ M© are 
expected to have e > 0.5. Again, the results are only 
mildly dependent on 7. Note that frequencies are com- 
puted in the reference frame of the source, the actual 
observed frequency has than to be appropriately red- 
shifted by a factor (1 -|- z) according to the redshift of 
the emitting system. This means that for high redshift 
sources, an ofeserwed frequency of 10"'' Hz corresponds to 
a higher intrinsic frequency. High redshift sources may 
then have milder eccentricities in the observable band, 
which may be relevant for LISA sources (whereas typical 
PTA sources are at z < 1). 

5.2. The impact of the chosen normalization 

We want to check at this point how the M — a normal- 
ization and the tuning of the loss cone refilling efficiency 
to Tinf impact on our results. Selected cases of alterna- 
tive models are presented in figure [7] The left hand plot 
is representative of LISA sources. We see that chang- 
ing a merely shifts the timescale of the binary evolution. 
For larger a, the system is more compact, the timescale 
for 3-body scattering is shorter, GW emission takes over 
later, and consequently the residual e at Jlisa is larger. 
However, as shown by the e(/fc) tracks, this is at most 
a factor of 2 effect. Changing the normalization p in 
the loss-cone refilling process, has instead a major im- 
pact on the evolutionary timescale and on the eccentric- 
ity evolution of the system, but still the residual eccen- 
tricity at Jlisa is basically unaffected when eo = 0.01, 
and it changes by at most a factor of 3 in the ep = 0.6 
case. In the right hand plot we consider instead the typ- 
ical PTA source. All the considerations made for the 
LISA case still hold, and in the relevant frequency range 
3 X 10~^Hz < / < 3 X 10~*Hz, the expected e changes by 
at most a factor of 2. We therefore consider our results 
quite robust irrespective of the assumed normalizations. 

5.3. Eccentricity distributions Jor selected MBHB 
population models 

As a final step, we quantify the eccentricity distribu- 
tion of GW sources in the relevant frequency range re- 
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Fig. 7. — Impact on the M — a assumption and on the p normalization on the evolution of the binary. Left plot: representative LISA 
source, with system parameter highlighted at the top. Right plot: typical PTA sources, with system parameters highlighted at the top. In 
each plot, in the left panels we considered three different normalizations of the stellar velocity dispersion: a = a (i.e., the value predicted 
by the M — a relation, solid lines), I.Sct (short-dashed lines), 0.76" (long— dashed lines); in the right panel we stick the efficiency of loss cone 
refilling to pj^f (solid lines), lOpinf (short-dashed lines), and O.lpinf (long-dashed lines). Top and middle panels represent the evolution of 
e and a against t, respectively; bottom panels represent e{fi^). In these latter panels, fhlSA and and the relevant frequency PTA range 
3 X IQ-^Hz < / < 3 X 10-*Hz are highlighted. 



suiting by applying our eccentricity evolution scheme to 
standard MBHB population models. 
For the LISA case we use two of the models utilized 



by the L ISA parameter estimation task force (Arun et. 
|al 1120091): in the first case seeds are light (M>1OOM0, 



model; Volonteri, Haardt fc Madau |2003), being 



the remnant o t the first Jr'UFiii star explosions (Madau 
2001|); in the second case, already quite heavy 



& Rees 



In this latter scenario, infact, sources are on average less 
massive and with low g, a condition that maximizes the 
eccentricity increase during the stellar scattering phase. 
If Co is already large (0.6 in our study case), then the 
observed eccentricity at /lis A is peaked at '-^ 0.1 for the 
BVR case and at ^ 0.4 in the VHM case. 

In exploring the consequences for PTA observations, 
we adopt the standard Tu-SA populati o n mo del em- 



(M>T(r Mq) seed BHs form by direct collapse of mas- 



(2009), where 



sive protogalactic discs (BVR model; Begelman, Volon- 
teri &: Rees [2006 1 . We ran 50 Monte Carlo realizations of 
each model, producing 50 catalogues of coalescing bina- 
ries over a period of 3 yrs. We then estimate the signal- 
to-noise ratio (SNR) of each binary in the LISA detector 
by assuming circular inspiral and computing the wave- 
form to the 2PN order. We then consider only those 
events resulting in an SNR> 8 in the detector, and we 
compute the expected eccentricity distribution at /lis A ■ 
Results are shown in figure [8] In the upper and in the 
middle panels we plot the contour plots of the differential 
distribution of GW sources as a function of Mi and q, 
cPN/dMidq, averaged over the 50 Monte Carlo realiza- 
tions, superposed to the contour plots for e\f^j^^{Mi,q). 
The two observed MBHB populations are extremely dif- 
ferent: in the VHM model, the bulk of sources have 
M ~ W^Mq with q - 0.1; while in the BVR model, 
most of the sources have M > 10'* Mq and g ss 1. The 
resulting eccentricities distributions are plotted in the 
lower panel, where we plot the probability density func- 
tion p(e) against e for the observed population. When 
eo = 0.01 (binaries approximately circular), eccentricity 
is expected to be < 10~^ in the BVR model, with a peak 
at about 2 x 10~^, but a broad eccentricity spectrum cov- 
ering the range 10^^ — 0.2 is expected in the VHM case. 



ployed by Sesana, Vecchio fc Volonteri 
merging galaxies are populated by MBl: 
the M — Afbuigc in the form given by Tundo et al. (12007 1 



s according to 



and accretion is triggered onto the more massive black 
ho le before the final coalescence. The r eader is referred 
to jSesana, Vecchio & Volonteri I (|2009|) for details. We 



ran 50 Monte Carlo realizations ot the model (assum- 
ing binaries in circular orbit) and we pick only the indi- 
vidually resolvable sources generating a timing residual 
larger than Ins. Again, the obtained differential distribu- 
tion of the individually resolvable sources, (PN/dMidq, 
averaged over the 50 Monte Carlo realizations, is su- 
perposed to the contour plots for e|pTA(-^i,9) in figure 
[9j The source distribution is strongly peaked around 
Ml — 10® M0 and q — 0.1, with a long tail extending to 
q — 10"'^. If binaries are approximately circular at the 
moment of pairing (eo = 0.01), then the expected p(e) is 
basically flat in the range [0.03,0.3], while for ep = (3.6, 
p(e) has a sharp peak in the range [0.5,0.7], highlighting 
the possible significant impact of eccentricity for PTA 
observations. 

6. DISCUSSION AND CONCLUSIONS 

We studied the semimajor axis and eccentricity evolu- 
tion of massive black hole binaries in stellar environments 
by coupling the results of numerical 3-body scattering ex- 
periments to an analytical framework describing the evo- 
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Fig. 8. — Evaluation of the eccentricity distribution of MBHBs 
observed by LISA. Top and middle panels, contour plots of the dif- 
ferential distribution of observable sources N / dM\dq (gray scale, 
contour normalization unnecessary for illustrative purposes), su- 
perimposed to the contour plots of e\f^^g^ (color scale; 7 = 1.5, 
eo = 0.01) in the {Mi,q) plane. Inset labels refer to the e con- 
tours. Top panel is for the VHM model, middle panel is for the 
BVR model. Bottom panel: probability density function p(e) cor- 
responding to the contour plot convolution; solid histograms are 
for eo = 0.01 and dashed histograms are for eo = 0.6. 
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Fig. 9. — Same as figure Is] but for PTA observations. In 
the top panel, the differential distribution of observable sources 
d^ N / dMidq is now superposed to e\f-p,^^, again assuming 7 = 1.5 
and eo = 0.01. The resulting p(e) is plotted in the bottom panel 
assuming eo = 0.01 (solid histogram) and eo = 0.6 (dashed his- 
togram) 

lution of the stellar distribution and the supply of stars 
to the binary loss cone. Our treatment takes into account 
the scattering of bound stars determining the erosion of 
the stellar cusp bound to the binary, and the subsequent 
scattering of unbound stars fed to the binary loss cone 
by relaxation processes. We do not address the nature 
of the relaxation processes leading to loss cone replenish- 
ment, but we treat the loss cone refilling efficiency as a 
parameter of the model. Eventually, GW emission takes 
over, leading to the final coalescence of the system. 

Our main finding is that 3-body scattering induces a 
significant increase in the MBHB eccentricity, that is not 
efficiently washed out by GW-induced circularization be- 
fore the system enter the LISA or the PTA bands. The 
eccentricity growth is in general larger for binaries with 
smaller mass ratios, and at the stellar scattering-GW 
transition can easily be higher than 0.9. Equal mass bi- 
naries in general experience a milder eccentricity growth 
when the initial eccentricity is close to zero. The ec- 
centricity growth is more prominent for systems charac- 
terized by smaller masses. Binaries with significant ini- 
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tial eccentricity eo > 0.3 end up in very eccentric orbits 
(e > 0.9) regardless on the other system parameters. 
The impact of the cusp slope can be significant, with 
shallower cusps leading to higher maximal values of e, as 
explained in SHM08. In general, the eccentricity growth 
is dominated by the bound scattering phase for binaries 
with q < 0.1 and by the unbound scattering phase for 
binaries with larger mass ratios. When compared to the 
sparse results of full N-body simulations found in liter- 
ature, the results of our models are in reasonable agree- 
ment with those of numerical studies. 

The implications for GW observations are relevant. 
When binaries are circular at the moment of pairing, 
their eccentricity when they enter the LISA band is in 
the range 10~^ — 0.2, and it is larger for low mass unequal 
binaries. If binaries are already eccentric at the moment 
of pairing, these figures shift to the range 10^'^ — 0.5, with 
lower mass binaries leading to higher eccentricities and 
only a mild dependence on the mass ratio. We empha- 
size once again that in our treatment, the total mass of 
the system sets the the typical scale of the problem; be- 
cause of this, the residual eccentricity in the LISA band 
is larger for lighter binaries. This is important because 
LISA will be mostly sensitive to low mass MBHBs in the 
range 10^ — 10^ M©. In the PTA windows, the implica- 
tions are even stronger. Initially circular systems end up 
with eccentricities in the range 10~^ — 0.8 at a frequency 
of 10~* Hz (relevant to PTA observations), and for sig- 
nificant initial eccentricity, binaries with Mi < 10^^ Mq 
always have e > 0.5 in the PTA band. The trend with the 
mass and the mass ratio are the same as for their LISA 
counterparts. All the results are basically independent 
on the cusp slope 7, and are only mildly dependent on 
the normalization of the stellar density distribution, and 
on the efficiency of the loss cone refilling. 

Once applied to standard MBHB population models, 
these results predict eccentricities in the range 10"'^ — 0.2 
(depending on the adopted seed formation model) for 
observable LISA sources, and a broad flat e distribution 
in the interval 0.03—0.3 for source individually resolvable 
by PTAs. High initial values of e, naturally lead to more 
eccentric systems. 

Our results are of particular interest for the GW com- 
munity, showing that a proper treatment of the eccen- 
tricity might be crucial in the challenge of GW detec- 
tion. Mock data challenge initiatives like the LISA mock 
data challenge, have so far implemented circular MBHBs 
only, and consequently, the ability of data analysis and 
parameter estimation algorithms has been proven only 
in this situation. The typical eccentricity values found 



in the LISA band (< 0.2 for systems in circular orbit 
at the moment of pairing) allow for a perturbative ap- 
proach to the problem of constructing trustwhorty post- 
Newtonian wavefor m, as the one recently employed by 
Yunes et al. (2009). In the light of the results presented 
here, further work in this direction would be extremely 
valuable. The addition of a non zero eccentricity would 
affect the waveform by adding significant amplitude mod- 
ulation and phase precession, which in turn would affect 
our detection and parameter estimation ability (work in 
this direc tion is ongoing and prelim inary results can be 



found in Porter fc Sesana (2010)). Also in the PTA 



source modelling held, the assumption of circular orbits 
ha s been widely used so far, wit h the notable exception 
of Enoki & NagashimaJ (20071. Further work on ec- 
centric source modelling is needed, in order to address 
properly how an eccentric population of MBHB would 
affect the overall level of the background, the statistics 
of individually resolvable sources, the detailed shape of 
the residuals and our ability of extracting signals and 
estimating source parameters. 

We finally stress that our model is oversimplified, rely- 
ing only on stellar dynamics without taking into account 
the possible impact of the presence of large amounts of 
gas surrounding the binary. Gas dynamics may be par- 
ticularly relevant to LISA sources, which are expected 
to be found in mergers of small ga laxies at high redshift 
( Sesana, Volonteri & Haardt |2007|) , where the mass con- 
tent of gala xies is likiely^to be dommated by gas (see, e.g.. 
Cole et al. |2000) . In this view, our model should provide 
a more trustful description of PTA sources wh ich con- 
sist instead of massive binaries at low redshift (Sesana, 



Vecchio &: Volonteri 2009), likely hosted by gas poorer 
galaxies. INonetheless, we should bear in mind that re- 
cent studies also found significant eccentri city increase 
in MBHBs driven by circumbinary dis ks (Armitage & 



Natarajan 2005 Cuadra et al. 2009). Moreover gas 



dynamics may not be efficient enough to drive the final 



coalescence of MBHB systems ( [Lodato et al. [2009 ) , and 
also for low mass sources at high redshift, stellar dynam- 
ics may provide a viable alternative path through the 
final coalescence. We hope that our exploratory study 
will stimulate further research on the subject, which is 
of critical importance for a comprehensive modelling of 
MBHB evolution and for their future observation in the 
upcoming gravitational radiation windows. 



I am grateful to E. Berti, M. Dotti and E. Porter for 
their comments and suggestions, and to Frank Ohme 
for the invaluable help in constructing the contour plots 
shown in the paper. 
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